Method and arrangement to detect and measure the phase of periodic bio-signals 

The invention relates to a method and an arrangement to reliably detect and measure in 
real time the phase of periodic physiological values or bio-signals. 

Methods are known according to the state of the art that use a staggered or sequentially- 
positioned analysis window across the temporal progression of the bio-signal to determine its 
phase. Methods or descriptive statistics based on Fourier transforms are applied to the signal 
section within the window. Thus, for example, periodically illuminating light marks of defined 
intensity and adequately-high frequency (above about 4 Hz) are used in order to monitor the 
functional capability of the visual system. An electro-encephalogram (EEG) is compiled for the 
function test, and the response to stimulus is analyzed with regard to the amplitudes and phase. 
The phase of the response to stimulus is one of the decisive diagnostic parameters in the 
functional diagnostic. 

A disadvantage of the conventional method is that the statistical unreliability of detection, 
or the inaccuracy of the measurement, is very high. This uncertainty and inaccuracy result from 
the signal theory as a result of, and in connection with, the length of the analysis window. The 
theory states that, as the length of the analysis window decreases, the statistical unreliability and 
thereby the inaccuracy increases, which is adequately known and has been adequately proved in 
practical signal analysis. The window must be enlarged in order to achieve a statistically 
improved result. It is known in the realm of physiology that the phase may alter relatively 
rapidly, and these alterations are also diagnostically relevant. With lengthier analyse windows 
the statistical unreliability of the measured result does not show the change and valuable 
information about phase changes is lost. 

It is the task of the invention to provide a method and an arrangement with which it is 
possible to detect and to measure the causal phase response in periodic bio-signals with better 
reliability and greater speed, and with simultaneous reduction in computing power, than when 
using conventional methods. 

This task is solved by the invention in that periodic bio-signals are determined 
corresponding to their physical and physiological source, in that a status observer is set up in 
parallel with the biological system under analysis, and in that a Kalman filter is used to evaluate 
the output values of the biological system and of the observer, and to determine the phase. 

In the method based on the invention, the phase of a periodic bio-signal is determined and 
used for functional diagnostic purposes. Thus, for example, a lengthened phase with respect to a 
healthy test subject may be an important clue to functional problems of the biological system 
being investigated. 

In the arrangement based on the invention, a status observer is arranged in parallel with 
the biological system being investigated that is imitated by a status model that, corresponding to 
system model, estimates the status value of phase based on a Kalman filter. 



Of advantage here are the facts that the estimation of the phase may occur continuously, 
and that no staggered or sequentially-applied analysis window is required. This makes analysis 
of the temporal phase alterations possible. In contrast to a relatively complicated theoretical 
background of this phase estimator, the practical implementation is simple. In comparison to the 
conventional method, it requires significantly reduced computing power, so that real-time phase 
estimation is possible. 

In the following, the invention will be described in greater detail using the theoretical 
derivation and an embodiment example. The pertinent Illustrations show: 

Figure 1 a flow chart of an observer concept; 

Figure 2 a system model of an arrangement based on the invention; 

Figure 3 a representation of the principle of the status observer to measure the 
phase in periodic bio-signals; 

Figure 4 a progression of an estimated phase for a harmonic of the frequency of 8 
Hz and Phase of 2 radians for the static Kalman factors 2 and 20; 

Figure 5 a progression of an estimated phase for a noise-affected harmonic of the 
frequency of 8 Hz and Phase of 2 radians with a SNR (signal-to-noise ratio) of 0 dB (lower) and 
the dynamic Kalman factor (upper); 

Figure 6 estimation of the signal's phase as in Figure 5 with static Kalman factors; 

Figure 7 results of the estimation of the phase (right) of real signals (left); 

Figure 8 an additive overlay of a harmonic of a frequency of 8 Hz with noise and 
intended de-tuning of the analysis frequency (above) and of the phase progression with 
acceleration (below). 

A biological system that produces a periodic bio-signal or responds to a periodic input 
signal is shown in Figure 1 as a status model of a "real system." The following status equations 
(1) and (2) describe this system (emboldened capital letters stand for matrices, and small letters 
stand for vectors): 



For further considerations, an additive signal model is assumed that sums a harmonic 
oscillation and normally-distributed noise: 



x(t) = A* x(t) + B • u(t); jc(0) = x(t 0 ) 



(1) 



y(t) = C*x(t). 



(2) 



y(t) = y*sm(co(t)) + r p (t) 



(3) 
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The goal is to construct a system model whose variable x(t) represents the phase (p(t) of 
the signal y(t) to be investigated. The phase cannot be measured directly since it is the argument 
of a trigonometric function. A supplemental construction is therefore required. One such 
construction is a status observer that is positioned in parallel with the system being investigated. 
The observer estimates the status variable by minimizing an error function that compares the 
outputs of the real system with those of the observer. In this manner, the status variable may be 
measured directly after successful error minimization. 

Figure 1 shows a flow chart of the observer concept. Since x(t) cannot be measured 
directly, x m (t) is estimated within the observer. The inner loop within the observer minimizes the 
error from y m (t) with respect to y(t) with the help of the correction matrix K. For the observer, 
the status equations (4) and (5) result: 

x M (t) = A.x M (t) + B.u(t) + K. [y(0 - y u (*)] , (4) 

y M {t) = Cx M (t). (5) 
From (4) and (5), we have: 

x M (t) = (A- KC) •x M (f) + B* «(/) + K . y{t) . (6) 

It is assumed that the systems possess different initial conditions. From this, we have the 
observation error: 

e(t) = x(t)-x M (t). (7) 

This observation error disappears iteratively with the help of the correction matrix K, 
resulting in, 

e(t) = 0 for t -> oo . (8) 

The dynamics and stability of the estimation may be described using the differential 
equation of the observer error (9): 

m=m- *„(*)■ (9) 

Simplification and additional intermediary steps lead to: 

e(t) = (A-K»Q»e(t). (10) 

Corresponding to the signal model (3), one must assume that the investigated signal is 
destroyed by noise. A Kalman filter is introduced to reduce the influence of noise. Taking the 
noise into account, the system is described by the following status equations: 

System status x(t) = A • x(t) + B» u{t) + M»r s (t) (11) 
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System output y(t) = C • x(t) + r p (t) ( 1 2) 

Observer x M (t) = A* x M (t) + B» u(t) + K(t) • \y(t) - C • x M (/)] , (13) 
whereby 

K(t) is the correction matrix that is to achieve the fact that e(t) = x{t) -x M (t)^0, 
e(t) = is the observer error, 
r s (t) is system noise, and 
r p (t) is process noise. 

In order to simplify the derivation, it is assumed that the noise components are wide-band 
Gaussian unit variance processes with known co-variances: 

cov rp {t x ,t 2 ) = E{r p ( tl ) • r T p (t 2 )}= R p (t } ) • 8{t x - 1 2 ) 

cov rj (f 1 ^ 2 ) = ^(r 1 )«^(^)}=^(^i)«^i-^) (14) 

the noise components are independent of one another, so that: 

cov Vs (t l ,t 2 ) = Q (15) 

For a constant estimation of x(t), the error performance must be minimized using the 
matrix K(t): 

E{e T {t) • e(t)} = E{ef(t) + e 2 2 (t) + -- + e 2 n (t)}= f(K(t)) = Min. , (16) 

Taking into account the stochastic relationships regarding the co-variances, a suitable 
correction matrix K(t) is derived corresponding to the Kalman filter: 

A-(0 = cov.(0»C r .J?; I (0. (17) 

The formula for the error co variance covs(t) may be derived from the Kalman filter: 

cov £ (0 = A • cov £ (t) + cov f (0 • A T - cov £ (f ) • C T • R' p \t) • cov £ (t) + M •R s (t)»M T . (18) 

Estimation of phase: 

The investigated signal is modeled based on (3) from the sum of a harmonic and the 



4 



v(0 = y„i (0 + r p (0 = y • sin(fi)/ + (pit)) + r p if), co = 2nf. (19) 
The phase results from the differential equation: 

<p(t) = -a*cp(t) + r s {t) a>0. (20) 

the system model shown in Figure 2 results from this. The phase cannot be measured 
directly. There is therefore parallel to the system an observer in which direct access to the 
estimated phase cp M (t) is possible. However, the non-linear component y n] (t) in (19) is 
unfavorable to the observer concept. A suitable linearization yi((p(t), t) is correspondingly 
required (5) in order create a linear relationship between the status variables cp M (t) and the output 
y M (t). Based on the Taylor linearization, the observer may be formulated as follows: 

<p M (0 =-a*(p M (0 + K(t) • CK0 - y M (0) • (21) 

^(0 = ^.(^(0.0- (22) 

The observer is modeled based on (21) and (22), as Figure 2 shows. Equation (22) is 
linked with Equation (5) in the result of the linearization at the working point cp B . The factor C 
that is used in (17) results from this in order to determine the correction factor K(t): 

K{t) = cov e it) • j> • cos(drt + (p B (0) • R'p (0 • (23) 

The phase to be determined is selected as the working point: 

P*(0 = M0. (24) 
and the resulting differential equation for the phase is: 

<Pm (0 = -a*<P M (0 + cov e (0 • y • cos(a* + (p M (0) • [y(0 - j> • sin(a* + <p M (/))]• ^' (0 . (25) 
In accordance with (25), the phase estimator may be modeled, as Figure 3 shows. 
In order to estimate the phase y(t), the error covariance must be calculated. From (18) 

results: 

cov e it) = -2a • cov e (/) - j> 2 • cos 2 (cot + (p M it)) • cov 2 (/) • R~ l (t) + R s (t) . (26) 

Equation (26) produces a simple solution if higher-frequency components are not taken 
into account in the error covariance. Based on (27), 

cos 2 icat + (p M it)) = 1 / 2 + cos(2atf + 2(p M it)) (27) 

equation 26 may be simplified to: 
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c o v e (t) = -2a • cov e (0 - 1 / 2 • y 2 • R~ l (t) • cov* (t) • +R S (t) . 



(28) 



Upon suitable selection of the parameter a in (28), high-frequency components are 
suppressed as the result of temporary integration, i.e., it possesses the properties of a low-pass 
filter. Taking the low-pass into account, Equation (25) may be simplified: 

<Pm (0 = -v<Pm (0 + cov e (0 • y • cos(<y/ + <p M (*)) • y(t) • R' p ] (0 , (29) 

Thus, the observer, shown in Figure 3, may be simplified. The system proposed in 
Equation (29) may particularly be used to estimate phase in harmonics within noise. 

Figure 4 shows the progression of the estimated phase for a harmonic at frequency 8 Hz 
and phase 2 radians for different Kalman factors. The Kalman factors are static, and equal either 
2 or 20. As may be taken from the graph, the estimation becomes slower as the Kalman factor is 
lower. Static Kalman factors must be used at the point where the point of the phase alteration in 
time is not known. 

Figure 5 shows the progression (lower part) of the estimated phase for a noise-influence 
harmonic of the frequency 8 Hz and phase 2 radians with a SNR of 0 dB and dynamic Kalman 
factor (upper part). If the point in time of the phase alteration is known, then the Kalman, factor 
may be so constructed that resultantly the alteration, and subsequently the variance of the 
estimation, may be reduced. 

In comparison, Figure 6 shows the phase estimation of the signal itself (in Figure 5) with 
static Kalman factors. 

Figure 7 shows the result of a phase estimation of real signals (right column of the 
graph). Sequences of light pulses with a repeat rate of 8 pulses per second are used to stimulate a 
visual system under investigation, for example, whereby the light stimulation is followed 
alternating with a rest pause (time of from 0 to 2 of the time progressions of the signals, upper 
left and lower left). The progression is shown of an EEG (electro-encephalogram, lift column of 
the graph) taken from two occipital positions after 16x stimulus-related averaging. Both time 
progressions show a clear jump in phase after the inception of light stimulation. 

Phase estimation becomes problematic for signals strongly influenced by noise. In 
general, it is true that the phase is more robust against interference than are the amplitudes, as is 
also known from the realm of Information Technology. However, in this fringe area, the 
question of the presence (or detection) of a causal phase must be determined, and only then may 
the phase be estimated. 

Figure 8 shows a harmonic or the frequency 8 Hz being additively overlaid with the noise 
beginning at t = 4 sec, whereby the SNR is -10 dB (upper curve in the graph). In this time 
interval, the amplitude of the harmonic cannot be established. If one applies the phase estimator 
with targeted de-tuning, here at a frequency of 7.8 Hz, or 0.2 Hz less than the frequency of the 
harmonic, then an increase of 1.2 radians/sec in the case of a causal phase (lower curve on the 
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graph). This increase may be used in combination with a discriminator directly to detect the 
signal. 



Reference symbol list 



a 


System parameters, selectable within the observer model 


A, B, C, K 


Matrices in the status model of a system 


9(t) 


Phase in the system model, value to be estimated 


<PM(t) 


Phase in the observer model, measurable value 


r P (t) 


Process noise 


r s (t) 


System noise 


u(t) 


Input variable of a system in the status model 


x(t) 


Status variable of a system in the status model 


x M (t) 


Status variable of the observer in the status model 


y(t) 


Output variable of a system in the status model 


y M (t) 


Output variable of the observer in the status model 


y.(9(t),t) 


Linearization operator for phase 
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